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Abstract 

A gap in phase-space, the loss cone (LC), is opened up by a supermassive black 
hole (MBH) as it disrupts or accretes stars in a galactic centre. If a star enters the LC then, 
depending on its properties, its interaction with the MBH will either generate a luminous 
electromagnetic flare or give rise to gravitational radiation, both of which are expected 
to have directly observable consequences. A thorough understanding of loss-cone refilling 
mechanisms is important for the prediction of astrophysical quantities, such as rates of tidal 
disrupting main-sequence stars, rates of capturing compact stellar remnants and timescales 
of merging binary MBHs. If a galaxy were isolated and perfectly spherical, the only refilling 
mechanism would be diffusion due to weak two-body encounters between stars. This would 
leave the LC always nearly empty. However, real galaxies are neither perfectly spherical nor 
isolated. In this thesis, we use A"-body simulations to investigate how noise from accreted 
satellites and other substructures in a galaxy's halo can affect the LC refilling rate. 

Any A r -body model suffers from Poisson noise which is similar to, but much stronger 
than, the two-body diffusion occurring in real galaxies. To lessen this spurious Poisson noise, 
we apply the idea of importance sampling to develop a new scheme for constructing A-body 
realizations of a galaxy model, in which interesting regions of phase-space are sampled by 
many low-mass particles. This scheme minimizes the mean-square formal errors of a given 
set of projections of the galaxy's phase-space distribution function. Tests show that the 
method works very well in practice, reducing the diffusion coefficients by a factor of ~ 100 
compared to the standard equal-mass models and reducing the spurious LC flux in isolated 
model galaxies to manageably low levels. 

We use multimass A-body models of galaxies with centrally-embedded MBHs to 
study the effects of satellite flybys on LC refilling rates. The total mass accreted by the MBH 
over the course of one flyby can be described, using a simple empirical fitting formula that 
depends on the satellite's mass and orbit. Published large-scale cosmological simulations 
yield predictions about the distribution of substructure in galaxy halos. We use results of 
these together with our empirical fitting formula to obtain an upper bound on substructure- 
driven LC refilling rates in real galaxies. We find that although the flux of stars into the 
initially emptied LC is enhanced, but the fuelling rate averaged over the entire subhalos 
is increased by only a factor 3 over the rate one expects from the Poisson noise due the 
discreteness of the stellar distribution. 
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Chapter 1 



Introduction 



In the high-redshift universe, massive black holes (MBHs) have been implicated as 
the powerhouses for quasars in the Active Galactic Nuclei (AGN) paradigm; in the nearby 
galaxies, the existence of dead quasar engines has also been supported by modern stellar- 
dynamical searches (e.g., Ferrarese & Ford 2005). Over three dozens of local detections 
have unveiled the demographics of MBH populations such as a tight correlation between 
MBH mass and galactic velocity dispersion (e.g., Gebhardt et al. 2000; Ferrarese <fc Merritt 
2000; Merritt & Ferrarese 2001), and a less tight correlation between MBH mass and the 
mass/luminosity of the bulge (e.g. Kormendy & Richstone 1995; Magorrian et al. 1998; 
Marconi & Hunt 2003). While efforts to build a larger and statistically significant sam- 
ple continue, it is time to understand the origin, evolution and cosmic relevance of these 
fascinating objects. 

The subject is large and my space is limited. In this thesis, I concentrate on 
developing and using iV-body models to learn more about the stellar dynamics in galactic 
centres, especially around the MBHs. I begin immediately below by reviewing the "journey" 
to the MBHs: from the AGN paradigm-mandated MBHs as power sources to the observation- 
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Chapter 1: Introduction 



supported MBHs. Then, I describe two astrophysical ingredients: electromagnetic flares 
by tidal disruption of main-sequence stars and gravitational wave signals by capture of 
compact stars. Next, I explain the emptiness of "loss cone" and the starvation of the MBH 
which motivate the work in thesis. Finally, I outline briefly iV-body modelling procedures 
developed in the subsequent chapters. 
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1.1 Supermassive black holes in galactic centres 

1.1.1 Energy arguments for MBHs in AGNs 

1.1.2 Loss of stars in the loss cone 

1.1.3 How to feed stars to a loss cone? 
Declaration 



Chapter 2 



Dynamical modelling 



Galaxies, as self-gravitating systems, can be idealized as configurations of point 
mass fully described by a phase-space distribution function (DF) f(x,v,t). The mass dis- 
tribution determines the gravitational field[ through Poisson's equation. Over their lifetime, 
galaxies are to a high degree collisionless; the DFs therefore satisfy the Collisionless Boltz- 
mann equation (Sec. I2.ip . With the help of Monte Carlo methods (Sec. I2.2.T]) . the evolution 
of the DF can be followed by iV-body integrations (Sec. I2.2."2l) . When encounters are taken 
into account, one writes the collision term in Master equation form (Sec. I2.3.T]) and expands 
in a Taylor series to derive the Fokker-Planck equation (Sec. 12.3^21) . 
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2.1 Equilibrium model 

2.1.1 DFs for collisionless stellar system 

2.1.2 Collisionless Boltzmann equation 

2.1.3 Preliminaries 
Units conversion and scalings 

2.2 A^body modelling 

2.2.1 Monte Carlo integration 

2.2.2 The CBE and iV-body simulations 

I. Initial conditions 

II. Poisson solvers 

III. Leap-frog integrator 
Summary 

2.3 Departure from equilibrium 

2.3.1 Collision terms 

2.3.2 Fokker-Planck equation 

2.3.3 Relaxation time 



Chapter 3 



Multi-mass schemes 

for collisionless iV-body simulations 



A general scheme for constructing Monte Carlo realizations of equilibrium, col- 
lisionless galaxy models with known distribution function (DF) /o is established. It uses 
importance sampling to find the sampling DF f s that minimizes the mean-square formal 
errors in a given set of projections of the DF /q. The result is a multi-mass iV-body realiza- 
tion of the galaxy model in which "interesting" regions of phase-space are densely populated 
by lots of low-mass particles, increasing the effective A there, and less interesting regions 
by fewer, higher-mass particles. 

The chapter is organized as follows. After recapping in the connection between 
A-body simulations and the CBE (Sec. 13.1.11) and presenting two possible weapons to fight 
against small A limitation (Sec. l3.1.2l) . we explain our multi-mass formulation in Sec. 13.21 In 
Sec. 13.31 we give two examples of using our scheme to suppress fluctuations in the monopole 
component of acceleration in spherical bulge models with or without a central MBH. We 
calculate formal estimates of the noise in A-body models constructed using equal-mass 
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scheme and SHQ95's method, and compare them to our own scheme (Sec. l3.3.2l) . In Sec. 13.3.31 
and 13.41 we test how well our realizations behave in practice when evolved using a real N- 
body code. For the bulge + MBH model in Sec. 13.41 much care has been taken to resolve the 
dynamics around the MBH. 
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3.1 Collisionless TV-body simulations 

3.1.1 iV-body simulations and the CBE 

3.1.2 How to fight against small iV limitations? 

3.2 Formulation 

3.2.1 Observables 

3.2.2 Optimal sampling scheme 

3.2.3 ICs for iV-body model 

3.3 Comparisons of a bulge model (no MBH) 

3.3.1 Formal errors 

I. The conventional equal-mass scheme 

II. Sigurdsson et al. 's multi-mass scheme 

III. Our scheme 
iV-body realizations 

3.3.2 How well is the acceleration field reproduced? 

3.3.3 How well are integrals of motion conserved? 

3.4 Simulations of a bulge + MBH model 

3.4.1 iV-body realizations 
II. Block leapfrog time-stepping 

3.4.2 Results: diffusion in J 2 



Chapter 4 



Loss cone refilling by flyby encounters 



In the last chapter, we have developed a general scheme for constructing multi-mass 
collisionless galaxy models. It has been successfully applied to set up a bulge + MBH model, 
whose long-time behavior has also been verified by full iV-body experiments. In this chapter, 
the MBH-embedded model is extended to include a loss cone (LC): a consumption sphere 
centred on the MBH. First of all, control tests have been directed to calibrate (artificial) 
noise-driven flux into the LC. Then starting from this near- equilibrium galaxy, satellites 
on flyby orbits are added as transient perturbers. For satellites with different mass and 
orbital parameters, I carry out extensive experiments to measure the mass of stars captured 
into the LC after a single fiyby encounter. Given a set of satellite parameters, an empirical 
relationship is found to predict the LC refilling mass. Finally, to understand the interactions 
between the perturber and stars, especially to disentangle effects of resonance coupling from 
velocity impulse, I run extra test-particle experiments of orbiting satellites. 
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4.1 Loss cone refilling by two-body relaxation 

What is the size of a loss cone? 

How to refill a loss cone in spherical galaxies? 

How to model a loss cone in iV-body simulations? 

4.2 Flyby encounters 

4.2.1 Tests: calibrating noise-driven LC refilling 

4.2.2 Experiments: measuring perturber-driven LC refilling 

I. Toy perturbers 

II. An example of flyby encounters 

4.2.3 Results: an empirical M\ c formula 
The validity of M\ c formula 

4.2.4 Discussion 

4.3 Where do LC stars come from? 

4.3.1 Angular-momentum transport by satellite torque 

4.3.2 Resonances from perturbation theory 

4.3.3 Are resonances important? 

Test-particle experiments of a periodic perturber 
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Rates of MBH feeding in bulge-halo 

systems 



In this chapter, I study the prospects for feeding the central MBH by orbiting 
halo substructures. I start by reviewing substructure properties observed in high-resolution 
cosmological iV-body simulations, then summarize information on the subhalo spatial dis- 
tribution, radial density profile, mass function and concentration parameters into a DF of 
the form f su b(£'i M). In Sec. 15.21 I use the subhalo DF together with the knowledge of how 
much mass can be fed to a MBH after a single encounter to estimate the perturber-driven 
LC refilling rate for a typical small galaxy. 
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Chapter 5: Rates of MBH feeding in bulge-halo systems 13 

5.1 Modelling a Milky- way sized halo and its substructures 

5.2 LC flux driven by noisy DM clumps 

5.3 Discussions 



Chapter 6 



Conclusions and future studies 



The motivation behind this thesis work was to improve our understanding of the 
internal dynamics at galactic centres, in particular, around super massive black holes. The 
best available tool on the dynamical study is undoubtedly full self-consistent iV-body simu- 
lations. We have (1) developed a multi-mass scheme for constructing collisionless iV-body 
models and (2) used the multi-mass bulge + MBH model to study loss cone dynamics. 
Below, I summarize results in earlier chapters and suggest some future work. 
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6.1 Constructing collisionless galaxy models 

In Chapter[2j I presented the basic dynamical equations that govern the evolution of 
collisionless systems. Beginning with the collisionless Boltzmann equation that is of essential 
importance for system (Sec, I2.1."2j) . I show explicitly how the CBE is related to the iV-body 
modelling (Sec. I2.2| ). When encounters are taken into account, one writes the collision term 
in Master equation form and expands in a Taylor series to derive the Fokker-Planck equation 
(Sec. El- 

Rather technical, Chapter [3] forms the backbone of the thesis. In it we saw how to 
design a sampling distribution function (DF) f s from some known DF /o (Sec. I3.2.2|) and 
how to generate a multi-mass collisionless model from it (Sec. I3.2.l3j) . In two applications, a 
bulge model and a bulge + MBH model, we aim at minimizing the shot noise in estimates 
of the acceleration field. Models constructed using our multi-mass scheme easily yield a 
factor ~ 100 reduction in the variance at the central acceleration field when compared 
to a traditional equal-mass model with the same number of particles. When evolving both 
models with an iV-body code, the diffusion coefficients in our model are reduced by a similar 
factor. Therefore, for certain types of problems, our scheme is a practical method for 
reducing the two-body relaxation effects, thereby bringing the iV-body simulations closer to 
the collisionless ideal. We note the following preparation for the applications of multi-mass 
modelling scheme: 

• For successful application, a system should be in a steady state, or close to one. 

• The DF /o should be quick and cheap to evaluate, either numerically or analytically. 
Finding /o for axisymmetric or triaxial galaxies is a longstanding and nontrivial prob- 
lem since one rarely has sufficient knowledge of the underlying potential's integrals 
of motion, but suitable flattened DFs do exist, including the standard axisymmetric 
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two-integral f(£,L z ) models and also rotating triaxial models such as those used in, 
e.g., Berczik et al. (2006) An alternative way of constructing flattened multi-mass 
realizations would be to apply Holley-Bockelmann et al. (2002) 's adiabatic sculpting 
scheme to a spherical A-body model constructed using our scheme. 

• The general multi-mass scheme uses importance sampling to find the tailored sampling 
DF f s that minimizes the sum of mean-square uncertainties in Qi (of the form eq. ??). 
As long as f s is smooth in integral space, Monte Carlo realizations of f s should work 
for any reasonably general collisionless A-body code. The utility of our multi-mass 
scheme, therefore, depends critically on the selection of the projection kernels Qi{w). 

The last point is new to this field. It is probably best addressed by experimenting with 
different sets of kernels, especially since it is easy to test the consequences of modifying 
them. Nevertheless, there are cases in which modest physical insight offers some guidance 
on choosing the Qi. Besides the loss cone problems addressed in Chapter 01 I here give 
another example for future work: 

^1. Sinking satellites Kazantzidis et al. (2004) demonstrate the significance of using 
equilibrium A-body realizations of satellite models when investigating the effect of 
tidal stripping of DM substructure halos (satellites) orbiting inside a more massive 
host potential. Besides the shape of the background potential and the amount of tidal 
heating, the mass-loss history is very sensitive to the detailed density profile of the 
satellite itself. One can therefore make one step further from equal-mass realizations 
by designing kernels to pick out orbits that pass through the tidal radius, while again 
maintaining an accurate estimate of the satellite's acceleration field. 
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6.2 Rates of loss cone refilling by MBHs 

Chapter H] describes applications of a multi-mass bulge + MBH model to study 
the loss cone (LC) dynamics. In Sec. 14.11 I introduced a loss cone consumption sphere 
to mimic tidal disruption events. The "empty" nature of the LC was then scrutinized by 
calibrating artificial LC refilling (Sec. I4.27T]) . Toy models in Sec. 14.2.21 taught us that only 
a fairly massive satellite (M s ) passing through the centre (b) with low speed (V) can give 
rise to a LC flux well above the threshold of numerical noise-driven LC refilling. Prom the 
reduced satellite parameter space, about 100 experiments have been carried out to measure 
the maximum net mass [M\ c ) captured into an initially empty LC. The collected data were 
summarized by a fitting Mi c -formula (??) in Sec. I4.2.3L 

• Stellar orbits become strongly perturbed during a satellite flyby. The LC-captured 
mass reaches a peak value soon after the perturber visits its pericentre. 

• M\ c is strongly correlated with the satellite parameters: scaling linearly with the satel- 
lite mass (M s ) and loss cone size (ri c ), but inversely with 6 2 ' and V 1 ' 3 . 

• As the satellite departs, stars settle down to a new equilibrium configuration within a 
few crossing times. Therefore, if perturbed by cumulative perturbers, one can simply 
accumulate each single event to get the total consumed mass. 

The rest of the Chapter [4] is devoted to understanding numerical results: where do LC 
stars come from? Linear perturbation theory claims that resonances are important if not 
most important in determining the response of a stellar system to any external perturbation 
(Sec. I4.3.2p . But for the flyby encounters of interest, could resonances be the culprit for 
redistributing angular-momentum in stars and eventually send stars to the LC? This was 
addressed in Sec. 14.3.31 where test-particle experiments were carried out to track resonant 
dynamics relevant to a periodic driver. I make the following summaries: 
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• Close encounters with the satellite effectively change angular-momentum of stars. The 
ones that lose enough angular-momentum can be effectively deflected into the loss cone 
and destroyed by falling into the MBH. 

• A flyby encounter does not repeat its pattern, which prevents stars in instantaneous 
resonance from building up a significant response. Therefore, resonance coupling be- 
tween stars and the external perturber is unimportant in terms of refilling the emptied 
loss cone. 

• The accelerated angular-momentum relaxation within a nearly Keplerian is also unim- 
portant here. 

In Chapter [H I introduced an outer halo as a reservoir to continuously inject orbiting subha- 
los into the inner bulge. One further step was made by describing the subhalo properties in 
terms of phase-space probability density, containing all the dynamical information (Sec l5.1|) . 
Finally in Sec. 15.21 I put things together to predict a LC refilling rate. The conclusion is: 

• The flux of stars into the loss cone is enhanced when the loss cone is initially emptied, 
but due to the scarcity of subhalo population near the galactic centre, the LC refilling 
rate averaged over the entire orbiting dark halo substructures is not strongly affected. 

The future 

For the galaxy model being used in the entire thesis, I have introduced several sim- 
plifications. At this point, I bring some notes about possible improvements and extensions 
of the current model. 

<s?2. Better galaxy models The galaxy model used in this thesis is not ideal: (a) the 
mass of the MBH is considerably larger than those observed in galactic centres and the 
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effective radius (r\ c ) of the LC is many orders of magnitude larger than the tidal radius 
in real galaxies. Unfortunately, it is still hard for simulations to determine accurately 
the flux of stars into a small MBH or a tiny LC. The only remedy is to increase the 
number of particles, the more the better, (b) The adopted Hernquist profile has 
a r^ 1 central cusp and no alternatives, such as core (shallower) or nucleus (steeper) 
profiles, have been explored here. It would be useful to extend the work in this thesis 
by building iV-body models with different mass profiles, MBH masses and LC sizes. 

^3. Evolution of binary supermassive black holes Larger elliptical galaxies and bulges 
grow through mergers. If more than one of the progenitors contains a MBH and 
the in-spiral time is less than a Hubble time, the MBHs will form a bound system. 
Begelman, Blandford & Rees (1980) showed that the evolution of binary MBHs in 
gas-poor galaxies can be divided into three phases: (a) As galaxies merge, the core 
undergoes violent relaxation. Via dynamical friction (Chandrasekhar 1943), the cap- 
tured MBH sinks towards the center of the common gravitational potential where they 
form a bound binary, (b) The binary continues to decay via gravitational slingshot 
interactions (Saslaw, Valtonen & Aarseth 1974): stars on orbits intersecting the bi- 
nary are ejected at much higher velocities comparable to the binary's orbital velocity, 
(c) Finally if the binary's separation decreases to the point at which gravitational 
radiation becomes the dominant dissipative force to carry away the last remaining 
energy and angular momentum, the binaries coalesce rapidly (Peters 1964). The long- 
term evolution of MBH binaries is always related to a long standing problem called 
the "final parsec problem" (Milosavljevic & Merritt 2003b): is super-elastic scattering 
off individual stars in the background efficient enough to transit a binary MBH from 
dynamical friction regime on the order of lpc scale, to the gravitational radiation dom- 
inated regime with a separation < 10~ 3 pc and finally merge within a Hubble time? 
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Uncertainties about the resolution have been a major impediment to predicting the 
frequency of MBH mergers in galactic nuclei, and hence to computing event rates for 
proposed gravitational wave interferometers like LISA. 1 As a natural continuation of 
the work in Chapter 03 I intend to further explore the use of Mi c -fitting-formula for 
the evolution of binary MBHs. 

■s?4. Response of galaxy to perturbations Our conclusion that resonances are unim- 
portant in redistributing angular-momentum within stars, is in disagreement with 
Vesperini & Weinberg (2000), who explore the effects of relatively weak encounters by 
computing the response of a spherical stellar system to the perturbation induced by 
a low-mass system during a flyby. They determine the contribution to the total self- 
consistent response from three sources: (a) the perturbation applied by the external 
perturber; (b) the reaction of the system to its own response to this perturbation; (c) 
the excitation of discrete damped modes of the primary system; and conclude that 
it is the resonances (point c) that lead to the excitation of patterns in the primary 
system. The main discrepancies probably stem from the galaxy model under study, 
since all their conclusions are based on a very strong assumption: discrete modes are 
weakly damped in King models with different concentrations. On the theoretical side, 
I intend to explore the response of galaxy to external perturbations by means of lin- 
ear perturbation theory, hopefully to find out whether weakly-damped modes exist in 
every galaxy. 



http://lisa.jpl.nasa.gov/ 
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